!********************
MODULE metric
!********************
USE ElliDef
USE NodeInfoDef
USE libinterp
IMPLICIT NONE
SAVE
PRIVATE

REAl (KIND=qPrec) :: deltax,deltay,deltaz,dmax,xmin,max,ymin,ymax,zmin,zmax
INTEGER nx,ny,nz

REAL (KIND=qPrec), DIMENSION (:,:), POINTER :: h,hx,hy,hxx,hyy

PUBLIC g,g_mask,setprob,DEFINE_TOPO



CONTAINS

SUBROUTINE setprob(Info)

TYPE (NodeInfo), TARGET :: Info

! set global variables in this module from values in Info

deltax=Info%dX(1);deltay=Info%dX(2);deltaz=Info%DX(3)
dmax=Info%dmax
nx=Info%mX(1);ny=Info%mX(2);nz=Info%mX(3)
xmin=Info%xLower(1);ymin=Info%xLower(2);zmin=Info%xLower(3)
IF(.NOT.ASSOCIATED(Info%topo)) ALLOCATE(Info%topo)
! set the local pointers to the topography in Info 
h=>Info%topo%h;hx=>Info%topo%hx;hy=>Info%topo%hy;
hxx=>Info%topo%hxx;hyy=>Info%topo%hyy;



END SUBROUTINE setprob

SUBROUTINE DEFINE_TOPO(Info,time)

TYPE (NodeInfo) :: Info
REAL (KIND=qPrec), OPTIONAL, INTENT (OUT) :: time
REAL (KIND=qPrec) :: dummy
REAl (KIND=qPRec), DIMENSION(:), ALLOCATABLE :: x,y,xg,yg
REAL (KIND=qPrec), DIMENSION (:), ALLOCATABLE :: w
REAL (KIND=qPrec), DIMENSION (:,:), ALLOCATABLE :: depth

INTEGER, DIMENSION (:), ALLOCATABLE :: iw
INTEGER :: lw, liw, intpol(2)=(/3,3/),ier,nxg,nyg,i,j

IF(PRESENT(time)) THEN
CALL CPU_TIME(dummy)
CALL CPU_TIME(dummy)
END IF


! interpolate topography
IF(.NOT.ASSOCIATED(Info%topo)) ALLOCATE(Info%topo)
ALLOCATE(Info%topo%h(0:Info%mX(1)+1,1),Info%topo%hx(0:Info%mX(1)+1,1))

h=>Info%topo%h;hx=>Info%topo%hx;hy=>Info%topo%hy;
hxx=>Info%topo%hxx;hyy=>Info%topo%hyy;


h=zero
! calculate derivatives
hx=zero
hx(1:Info%mX(1),1)=(h(2:Info%mX(1)+1,1)-h(0:Info%mX(1)-2,1))/(2*Info%dX(1))
! zero the boundaries, we do not like them nonzero
hx(0,:)=(h(1,:)-h(0,:))/(Info%dX(1));
hx(Info%mX(1)+1,:)=(h(Info%mX(1)+1,:)-h(Info%mX(1),:))/(Info%dX(1))
PRINT*,minval(hx),'<hx<',maxval(hx)

 


IF(PRESENT(time)) THEN
   CALL CPU_TIME(time)
time=time-dummy
END IF


END SUBROUTINE DEFINE_TOPO


ELEMENTAL REAL (KIND=qPrec) FUNCTION g(i,j,k,l,p)
! returns \sqrt(g)glp at location i,j,k
! NOTE THAT WHEN l==p IT IS RETURNED AT THE CORRESPONDING
! CELL EDGE
INTEGER, INTENT (IN) :: i,j,k,l,p

REAL (KIND=qPRec) :: xi,eta,zeta




SELECT CASE( (l-1)*2+p)
CASE (-2) !\sqrt(g)
g=(1.-h(i,1))

CASE (1) !g11
xi=(i-1)*deltax+xmin
eta=(j-.5)*deltay+ymin
g=(1.-(h(i,1)+h(i-1,1))/(2))

CASE (4) !g22
xi=(i-.5)*deltax+xmin
eta=(j-1)*deltay+ymin


g=(1.+(eta-1)**2*(hx(i,1)**2))/(1.-h(i,1))

CASE (2,3) ! g12

xi=(i-.5)*deltax+xmin
eta=(j-.5)*deltay+ymin


g=(eta-1)*hx(i,1)

END SELECT
END FUNCTION g



LOGICAL FUNCTION g_mask(i,j)

INTEGER, INTENT (IN) :: i,j
SELECT CASE ((i-1)*2+j)
CASE(1) ! 11
g_mask=.true.   
CASE(4) ! 22
g_mask=.true.   
CASE(2,3) ! 12 or 21
g_mask=.true.   
CASE DEFAULT
   PRINT*,'g_mask not defined for',i,j
   STOP
END SELECT
END FUNCTION g_mask
END MODULE metric
